/*
 ***************************************************************************
 *
 * \file   hoffman2.cpp
 * Laplace equation in 2D
 *
 * \par Description
 * This code solves the Laplace equation in two dimensions: 
 * \f$ \nabla^2 T = 0\f$. 
 * \par
 * This equation is solved in the rectangle
 * interval \f$ [0,1] \times [0,2] \f$. The boundary conditions are ...
 * \par
 *  To compile and run this example type the next commands:
 * \par
   \verbatim
   % make
   % ./hoffman2 \endverbatim
 * 
 ***************************************************************************
 *  \author    Luis M. de la Cruz Salas [ Fri Mar 21 12:26:33 GMT 2008 ]
 ***************************************************************************
 */

#include <iostream>
#include <cmath>
#include "Meshes/StructuredMesh.hpp"
#include "Meshes/Uniform.hpp"
#include "Equations/ScalarEquation.hpp"
#include "Schemes/LaplaceScheme.hpp"
#include "Solvers/TDMA.hpp"
#include "Utils/num_utils.hpp"
#include "Utils/inout.hpp" 
#include "Utils/GNUplot.hpp"
using namespace Tuna;

// Scalar field in two dimension
typedef TunaArray<double, 2>::huge ScalarField2D;

void exactSolution(ScalarField2D &T, double hx, double hy, double A,
		   double dx, double dy)
{
  int NI = T.ubound(firstDim);
  int NJ = T.ubound(secondDim);
  
  for(int i = 1; i < NI; ++i)
    for(int j = 1; j < NJ; ++j) 
      T(i,j) = A * sinh( PI * (hy - j * dy) / hx ) * sin( PI * i * dx / hx) / sinh( PI * hy / hx);     
}

int main() 
{
  double length_x, length_y, Amp;
  int num_nodes_x, num_nodes_y;
  double tolerance;    
  int tdma_iter;
  timer time;
  
// Reading data from the input file
  std::ifstream input_cin ("input_hoffman"); 
  input_cin >> length_x;
  input_cin >> length_y;
  input_cin >> num_nodes_x; 
  input_cin >> num_nodes_y; 
  input_cin >> Amp; 
  input_cin >> tdma_iter;
  input_cin >> tolerance;
  input_cin.close();

  std::cout.precision(10);
  std::cout.setf(ios_base::fixed,ios_base::floatfield);

  StructuredMesh<Uniform<double, 2> > mesh(length_x, num_nodes_x,
				           length_y, num_nodes_y);
  double dx = mesh.getDelta(X);
  double dy = mesh.getDelta(Y);
  mesh.print();

  ScalarField2D T (mesh.getExtentVolumes());
  ScalarField2D Ta(mesh.getExtentVolumes());
  ScalarField2D Tn(mesh.getExtentNodes());
  
  Range all = Range::all();

  firstIndex i;
  Range I(Tn.lbound(firstDim), Tn.ubound(firstDim));
  Tn(I, Tn.lbound(secondDim)) = 100 * sin( PI * dx * i/ length_x);
  InOut::writeToFile(Tn, 0, "./DataHoffman2/temp.", dx, dy);
  NumUtils::interpolateToVolumes(T, Tn);
  
  exactSolution(Tn, length_x, length_y, Amp, dx, dy);
  InOut::writeToFile(Tn, 0, "./DataHoffman2/exac.", dx, dy);
  NumUtils::interpolateToVolumes(Ta, Tn);
 
#ifdef WITH_GNUPLOT
    GNUplot plotter;
    plotter("set grid");
    plotter("set contour");
    plotter("set view 60,45,1,1");
    plotter("set ticslevel 1.5");
    plotter("set cntrparam level discrete 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95");
    plotter("splot [0:1][0:2][0:100] \"./DataHoffman2/exac.00000\" title \"Exact solution\" w l");
    std::cout << "\n >---> Press <1> and then <enter> to continue \n";
    int pause; std::cin >> pause;
    plotter("splot [0:1][0:2][0:100] \"./DataHoffman2/temp.00000\" title \"Initial conditions\"w l");
    std::cout << "\n >---> Press <1> and then <enter> to continue \n";
    std::cin >> pause;
#endif

    SparseMatrix< Diagonal<double, 2> > A(num_nodes_x, num_nodes_y); 

    ScalarField2D             b(num_nodes_x, num_nodes_y);
  
    ScalarEquation< LaplaceScheme<double, 2> > laplace(T, A, b,
						       mesh.getDeltas());
    laplace.setDirichlet(TOP_WALL);
    laplace.setDirichlet(BOTTOM_WALL);
    laplace.setDirichlet(LEFT_WALL);
    laplace.setDirichlet(RIGHT_WALL);
    laplace.setGamma(1.0);
    laplace.print();

    time.tic();        

    laplace.calcCoefficients();
    int iter = Solver::TDMA2DY(laplace, tolerance, tdma_iter);

    std::cout << "\n Elapsed time for solving Ax = b : " 
	      << time.toc() << " seconds \n\n";   
    
    std::cout << "\n Iterations = " << iter
	      << "\t Residual = " << laplace.getResidual();

    laplace.update();  
    NumUtils::interpolateToNodes(Tn, T);
    InOut::writeToFile(Tn, 1, "./DataHoffman2/temp.", dx, dy);

    cout << "\t Error = " << NumUtils::calcErrorL1( T, Ta ) << endl;    
       
    for(int i = 0; i <= T.ubound(firstDim); ++i)
      for(int j = 0; j <= T.ubound(secondDim); ++j) 
	T(i,j) = abs(T(i,j) - Ta(i,j));




    NumUtils::interpolateToNodes(Tn, T);           
    InOut::writeToFile(Tn, 0, "./DataHoffman2/erro.", dx, dy);

#ifdef WITH_GNUPLOT
    plotter("splot [0:1][0:2][0:100] \"./DataHoffman2/temp.00001\" title \"Numerical Solution\" w l");
    std::cout << "\n >---> Press <1> and then <enter> to continue \n";
    std::cin >> pause;
    plotter("splot [0:1][0:2][0:100] \"./DataHoffman2/erro.00000\" title \"Error\" w l");
    std::cout << "\n >---> Press <1> and then <enter> to continue \n";
    std::cin >> pause;
#endif
    
    return 0;
}




